PROGRAM Hess_Smith                                              ! my English is not good , existing some grammatical errors
    USE CPUTime_MOD
    IMPLICIT NONE
    REAL(8) :: A1, B1                                ! for calculating the volume of the ellipse
    Integer(4) :: M                                    ! famous constant
    REAL(8), ALLOCATABLE :: Panel(:, :, :)                         ! panel points location
    REAL(8) :: Volume, Pi = 3.1415926
    REAL(8) :: T, I, J, K
    REAL(8), ALLOCATABLE :: Panel_Normal(:, :)
    REAL(8), ALLOCATABLE :: Panel_Centroid(:, :)
    REAL(8), ALLOCATABLE :: Loc_Coord(:, :, :)
    REAL(8), ALLOCATABLE :: Loc_Point(:, :, :)
    REAL(8), ALLOCATABLE :: Panel_Area(:)
    REAL(8), ALLOCATABLE :: Panel_N(:, :)
    REAL(8), ALLOCATABLE :: A(:, :)
    REAL(8), ALLOCATABLE :: B(:, :)
    REAL(8), ALLOCATABLE :: U(:, :)                               ! distribution source intensity
    REAL(8), ALLOCATABLE :: S(:, :)
    REAL(8), ALLOCATABLE :: Phi(:, :)
    REAL(8) :: Add_Mass(6, 6)

    CALL MESH(A1, B1, M)
    ! WRITE(*,*) "A,B,M"
    ! READ (*,*) A1,B1,M
    OPEN (10, File='Input.Dat', STATUS='Old', ACTION='READ')      ! input the panel information ; 12 vertical lines ; include the location of the penel points (4)
    ALLOCATE (Panel(M, 4, 3))
    READ (10, *) (((Panel(I, J, K), K=1, 3), J=1, 4), I=1, M)       ! input all the panel nodes
    CLOSE (10)

    Call GetStartTime()
    Volume = Pi*4/3*A1*B1**2

    ALLOCATE (Panel_Normal(M, 3), Panel_Centroid(M, 3), Loc_Point(M, 4, 3), Loc_Coord(M, 3, 3), Panel_Area(M))
    CALL Local(M, Panel, Panel_Normal, Panel_Centroid, Loc_Point, Loc_Coord, Panel_Area)
    ALLOCATE (A(M, M))

    DO J = 1, M
    DO I = 1, M
        IF (I == J) THEN
            A(I, J) = 2*Pi
        Else
            CALL Source_Distribution(Panel_Centroid(J,1:3),Panel_Centroid(I,1:3),Loc_Point(J,1:4,1:3),Loc_Coord(I,3,1:3),Loc_Coord(J,1:3,1:3),A(I,J))
        END IF
    END DO
    END DO

    ALLOCATE (B(M, 6), Panel_N(M, 6))

    DO I = 1, M
        CALL Source_B(Loc_Coord(I, 3, 1:3), Panel_Centroid(I, 1:3), B(I, 1:6), Panel_N(I, 1:6))
    END DO

    ALLOCATE (U(M, 6))

    CALL Brinv(A, M)
    U = MATMUL(A, B)
    ALLOCATE (S(M, M))

    DO J = 1, M
        DO I = 1, M
            CALL Potential(Panel_Centroid(J, 1:3), Panel_Centroid(I, 1:3), Loc_Point(J, 1:4, 1:3), Loc_Coord(J, 1:3, 1:3), S(I, J))
        END DO
    END DO

    ALLOCATE (Phi(M, 6))                                        ! Phi is the velocity potential

    DO K = 1, 6
        DO I = 1, M
            T = 0
            DO J = 1, M
                T = T + S(I, J)*U(J, K)
            END DO
            Phi(I, K) = T                                        ! Phi[M,6] means the 6-dimensional potentials in the ith panel
        END DO
    END DO

    DO I = 1, 6
        DO J = 1, 6
            T = 0
            DO K = 1, M
                T = T + Panel_Area(K)*Panel_N(K, J)*Phi(K, I)
            END DO
            Add_Mass(J, I) = T/Volume
        END DO
    END DO

    Call GetEndTime()
    OPEN (30, File='Output.Dat')
    WRITE (30, *) "A,B,M,Volume"
    WRITE (30, *) A1, B1, M, Volume
    WRITE (30, *) '# [ Add Mass Coefficient Matrix ] '
    WRITE (30, "(6(1X,ES9.2))") ((Add_Mass(J, I), I=1, 6), J=1, 6)
    WRITE (30, *) "I,Loc_Coord(I,3,:),Panel_Area(I),Phi(I,1:3)"
    WRITE (30, 101) ((I, Loc_Coord(I, 3, :), Panel_Area(I), Phi(I, 1:3)), I=1, M)
101 FORMAT(8(1X, ES9.2))

    WRITE (*, *) "Calculation over, press any key to quit..."
    READ (*, *)

END PROGRAM

SUBROUTINE Local(M, Panel, Panel_Normal, Panel_Centroid, Loc_Point, Loc_Coord, Panel_Area)
    IMPLICIT NONE
    Integer(4) :: M
    REAL(8) :: Panel(M, 4, 3)
    Integer(4) :: I, J, K
    REAL(8) :: Panel_Normal(M, 3)
    REAL(8) :: Panel_Centroid(M, 3)
    REAL(8) :: Loc_Point(M, 4, 3)
    REAL(8) :: Loc_Coord(M, 3, 3)
    REAL(8) :: Area_Point(M, 4, 3)
    REAL(8) :: Panel_Area(M)
    REAL(8) :: L, L1, L2, L3, L4, P1, P2, P
    REAL(8) :: Tri(M, 2, 3)
    REAL(8) :: Tri_Area(M, 2)
    REAL(8) :: N1, N2, N3, M1, M2, M3

    DO I = 1, M
        Tri(I, 1, 1) = (Panel(I, 1, 1) + Panel(I, 2, 1) + Panel(I, 3, 1))/3                                    ! convert one rectangle to two triangles
        Tri(I, 1, 2) = (Panel(I, 1, 2) + Panel(I, 2, 2) + Panel(I, 3, 2))/3
        Tri(I, 1, 3) = (Panel(I, 1, 3) + Panel(I, 2, 3) + Panel(I, 3, 3))/3
        Tri(I, 2, 1) = (Panel(I, 1, 1) + Panel(I, 4, 1) + Panel(I, 3, 1))/3                                    ! maybe triangle becomes a line since panel.node.3 = panel.node.4
        Tri(I, 2, 2) = (Panel(I, 1, 2) + Panel(I, 4, 2) + Panel(I, 3, 2))/3
        Tri(I, 2, 3) = (Panel(I, 1, 3) + Panel(I, 4, 3) + Panel(I, 3, 3))/3
    END DO

    DO I = 1, M
L1 = SQRT((Panel(I, 1, 1) - Panel(I, 2, 1))**2 + (Panel(I, 1, 2) - Panel(I, 2, 2))**2 + (Panel(I, 1, 3) - Panel(I, 2, 3))**2)        ! panel nodes' distances (4)
       L2 = SQRT((Panel(I, 2, 1) - Panel(I, 3, 1))**2 + (Panel(I, 2, 2) - Panel(I, 3, 2))**2 + (Panel(I, 2, 3) - Panel(I, 3, 3))**2)
       L3 = SQRT((Panel(I, 3, 1) - Panel(I, 4, 1))**2 + (Panel(I, 3, 2) - Panel(I, 4, 2))**2 + (Panel(I, 3, 3) - Panel(I, 4, 3))**2)
       L4 = SQRT((Panel(I, 4, 1) - Panel(I, 1, 1))**2 + (Panel(I, 4, 2) - Panel(I, 1, 2))**2 + (Panel(I, 4, 3) - Panel(I, 1, 3))**2)
L = SQRT((Panel(I, 3, 1) - Panel(I, 1, 1))**2 + (Panel(I, 3, 2) - Panel(I, 1, 2))**2 + (Panel(I, 3, 3) - Panel(I, 1, 3))**2)        ! panel's diagonal length
        P1 = (L1 + L2 + L)/2
        P2 = (L3 + L4 + L)/2
        Tri_Area(I, 1) = SQRT(P1*(P1 - L1)*(P1 - L2)*(P1 - L))                                            ! Heron's formula , calculate the area of a triangle : S = SQRT( p*(p-a)*(p-b)*(p-c) )
        Tri_Area(I, 2) = SQRT(P2*(P2 - L3)*(P2 - L4)*(P2 - L))
        Panel_Area(I) = Tri_Area(I, 1) + Tri_Area(I, 2)                                                ! merge two triangles to one rectangle (area)
    END DO

    DO I = 1, M                                                                                     ! calculate the (x,y,z) of panels' centroids , J=1=x , J=2=y , J=3=z
        DO J = 1, 3                                                                                 ! centroid coordinate  calculate formula
            Panel_Centroid(I, J) = (Tri(I, 1, J)*Tri_Area(I, 1) + Tri(I, 2, J)*Tri_Area(I, 2))/Panel_Area(I)
        END DO
    END DO

    DO I = 1, M
        N1 = Panel(I, 3, 1) - Panel(I, 1, 1)                                                             ! a panel's 3 dimensional distances two diagonal lengths
        N2 = Panel(I, 3, 2) - Panel(I, 1, 2)
        N3 = Panel(I, 3, 3) - Panel(I, 1, 3)
        M1 = Panel(I, 2, 1) - Panel(I, 4, 1)
        M2 = Panel(I, 2, 2) - Panel(I, 4, 2)
        M3 = Panel(I, 2, 3) - Panel(I, 4, 3)
        Panel_Normal(I, 1) = N2*M3 - N3*M2                                                            ! a panel's normal vector
        Panel_Normal(I, 2) = N3*M1 - N1*M3
        Panel_Normal(I, 3) = N1*M2 - N2*M1
        P = SQRT(Panel_Normal(I, 1)**2 + Panel_Normal(I, 2)**2 + Panel_Normal(I, 3)**2)
        Panel_Normal(I, 1) = Panel_Normal(I, 1)/P                                                    ! calculate the unit normal vector
        Panel_Normal(I, 2) = Panel_Normal(I, 2)/P
        Panel_Normal(I, 3) = Panel_Normal(I, 3)/P
    END DO

    DO I = 1, M
        Loc_Coord(I, 3, 1) = Panel_Normal(I, 1)                                                       ! it is a unit vector now
        Loc_Coord(I, 3, 2) = Panel_Normal(I, 2)
        Loc_Coord(I, 3, 3) = Panel_Normal(I, 3)

        DO J = 1, 4
            K = DOT_PRODUCT(Loc_Coord(I, 3, :), Panel(I, J, :) - Panel_Centroid(I, :))                   ! form a surface element throught the Panel_Centroid whose normal vector is Panel_Normal
            Loc_Point(I, J, :) = Panel(I, J, :) - K*Loc_Coord(I, 3, :)
        END DO
        Loc_Coord(I, 1, :) = Loc_Point(I, 1, :) - Loc_Point(I, 3, :)                                       ! ? form a unit vector in the surface , applicable to rectangle and triangle elements
        P = SQRT(Loc_Coord(I, 1, 1)**2 + Loc_Coord(I, 1, 2)**2 + Loc_Coord(I, 1, 3)**2)
        Loc_Coord(I, 1, :) = Loc_Coord(I, 1, :)/P                                                      ! unitization
        Loc_Coord(I, 2, 1) = Loc_Coord(I, 1, 2)*Loc_Coord(I, 3, 3) - Loc_Coord(I, 1, 3)*Loc_Coord(I, 3, 2)     ! cross product
        Loc_Coord(I, 2, 2) = Loc_Coord(I, 1, 3)*Loc_Coord(I, 3, 1) - Loc_Coord(I, 1, 1)*Loc_Coord(I, 3, 3)
        Loc_Coord(I, 2, 3) = Loc_Coord(I, 1, 1)*Loc_Coord(I, 3, 2) - Loc_Coord(I, 1, 2)*Loc_Coord(I, 3, 1)

        DO J = 1, 4
            N1 = Panel(I, J, 1) - Panel_Centroid(I, 1)                                         ! the 3-dimensional distances between the centroids of the ith and jth Panel element
            N2 = Panel(I, J, 2) - Panel_Centroid(I, 2)
            N3 = Panel(I, J, 3) - Panel_Centroid(I, 3)
            Loc_Point(I, J, 1) = N1*Loc_Coord(I, 1, 1) + N2*Loc_Coord(I, 1, 2) + N3*Loc_Coord(I, 1, 3)                     ! x
            Loc_Point(I, J, 2) = N1*Loc_Coord(I, 2, 1) + N2*Loc_Coord(I, 2, 2) + N3*Loc_Coord(I, 2, 3)                     ! y
            Loc_Point(I, J, 3) = N1*Loc_Coord(I, 3, 1) + N2*Loc_Coord(I, 3, 2) + N3*Loc_Coord(I, 3, 3)                     ! z
        END DO

    END DO

    RETURN
END SUBROUTINE Local

SUBROUTINE Source_Distribution(Panel_Centroid_J, Panel_Centroid_I, Loc_Point, Loc_Coord_IZ, Loc_Coord, A)
    ! function :
    ! to calculate influence coefficients , which mean the influence on the ith control point from the distributed source on the jth unit
    IMPLICIT NONE
    REAL(8) :: Panel_Centroid_I(3), Panel_Centroid_J(3)
    REAL(8) :: Loc_Point(4, 3), Loc_Coord_IZ(3)
    REAL(8) :: Loc_Coord(3, 3)
    Integer(4) :: I, J, K
    REAL(8) :: S, A, B, A1, B1
    REAL(8) :: Loc_P(3)
    REAL(8) :: Qq(5, 3)
    REAL(8) :: Sx, Sy, Sz, L, R1, R2, C1, C2, H1, H2, Mm
    REAL(8) :: N1, N2, N3, M1, M2, M3

    DO I = 1, 5
        DO J = 1, 3
            IF (I <= 4) THEN
                Qq(I, J) = Loc_Point(I, J)                                                 ! deliver the nodes location to Qq in the local coordinate
            Else
                Qq(5, J) = Loc_Point(1, J)
            END IF
        END DO
    END DO

    N1 = Panel_Centroid_I(1) - Panel_Centroid_J(1)                                         ! the 3-dimensional distances between the centroids of the ith and jth Panel element
    N2 = Panel_Centroid_I(2) - Panel_Centroid_J(2)
    N3 = Panel_Centroid_I(3) - Panel_Centroid_J(3)
    Loc_P(1) = N1*Loc_Coord(1, 1) + N2*Loc_Coord(1, 2) + N3*Loc_Coord(1, 3)                     ! the location
    Loc_P(2) = N1*Loc_Coord(2, 1) + N2*Loc_Coord(2, 2) + N3*Loc_Coord(2, 3)
    Loc_P(3) = N1*Loc_Coord(3, 1) + N2*Loc_Coord(3, 2) + N3*Loc_Coord(3, 3)

    A = 0
    B = 0
    S = 0

    DO I = 1, 4
        L = SQRT((Qq(I + 1, 1) - Qq(I, 1))**2 + (Qq(I + 1, 2) - Qq(I, 2))**2 + (Qq(I + 1, 3) - Qq(I, 3))**2)            ! distances of nearby nodes
        R1 = SQRT((Loc_P(1) - Qq(I, 1))**2 + (Loc_P(2) - Qq(I, 2))**2 + (Loc_P(3) - Qq(I, 3))**2)               ! Loc_P is the node not in the local panel
        R2 = SQRT((Loc_P(1) - Qq(I + 1, 1))**2 + (Loc_P(2) - Qq(I + 1, 2))**2 + (Loc_P(3) - Qq(I + 1, 3))**2)
        IF (L == 0) THEN
            A1 = 0
            B1 = 0
        Else
            A1 = -(Qq(I + 1, 2) - Qq(I, 2))/L*Log((R1 + R2 + L)/(R1 + R2 - L))
            B1 = (Qq(I + 1, 1) - Qq(I, 1))/L*Log((R1 + R2 + L)/(R1 + R2 - L))
        END IF
        A = A + A1
        B = B + B1
    END DO

    Sx = A
    Sy = B

    DO I = 1, 4
        C1 = (Qq(I, 1) - Loc_P(1))**2 + Loc_P(3)**2
        C2 = (Qq(I + 1, 1) - Loc_P(1))**2 + Loc_P(3)**2
        H1 = (Qq(I, 1) - Loc_P(1))*(Qq(I, 2) - Loc_P(2))
        H2 = (Qq(I + 1, 1) - Loc_P(1))*(Qq(I + 1, 2) - Loc_P(2))
        IF (Qq(I, 1) == Qq(I + 1, 1)) THEN
            Mm = 0
        Else
            Mm = (Qq(I + 1, 2) - Qq(I, 2))/(Qq(I + 1, 1) - Qq(I, 1))
        END IF
        R1 = SQRT(C1 + (Qq(I, 2) - Loc_P(2))**2)
        R2 = SQRT(C2 + (Qq(I + 1, 2) - Loc_P(2))**2)
        S = S + Atan((Mm*C1 - H1)/(Loc_P(3)*R1)) - Atan((Mm*C2 - H2)/(Loc_P(3)*R2))
    END DO

    Sz = S
    A = Loc_Coord_IZ(1)*Sx + Loc_Coord_IZ(2)*Sy + Loc_Coord_IZ(3)*Sz

    RETURN
END SUBROUTINE Source_Distribution

SUBROUTINE Source_B(Panel_Normal, Panel_Centroid, B, N)
    IMPLICIT NONE
    REAL(8) :: Panel_Normal(3)
    REAL(8) :: Panel_Centroid(3)
    REAL(8) :: P
    REAL(8) :: B(6), N(6)
    Integer(4) :: I, J, K
    REAL(8) :: R(3)

    P = SQRT(Panel_Centroid(1)**2 + Panel_Centroid(2)**2 + Panel_Centroid(3)**2)
    R(1) = Panel_Centroid(1)/P
    R(2) = Panel_Centroid(2)/P      ! unitization
    R(3) = Panel_Centroid(3)/P
    DO I = 1, 3
        N(I) = Panel_Normal(I)
    END DO

    N(4) = R(2)*N(3) - R(3)*N(2)                !  vector R X vecter N
    N(5) = R(3)*N(1) - R(1)*N(3)
    N(6) = R(1)*N(2) - R(2)*N(1)

    DO I = 1, 6
        B(I) = N(I)                           ! this place, V is unit number.
    END DO

    RETURN
END SUBROUTINE Source_B

SUBROUTINE Potential(Panel_Centroid_J, Panel_Centroid_I, Loc_Point, Loc_Coord, S)
    IMPLICIT NONE
    REAL(8) :: Panel_Centroid_I(3), Panel_Centroid_J(3)
    REAL(8) :: Loc_Point(4, 3)
    REAL(8) :: Loc_Coord(3, 3)
    Integer(4) :: I, J
    REAL(8) :: S, S1
    REAL(8) :: Loc_P(3)
    REAL(8) :: Qq(5, 3)
    REAL(8) :: Sz
    REAL(8) :: L, R1, R2, C1, C2, H1, H2, Mm
    REAL(8) :: N1, N2, N3

    DO I = 1, 5
        DO J = 1, 3
            IF (I <= 4) THEN
                Qq(I, J) = Loc_Point(I, J)                                                 ! deliver the nodes location to Qq in the local coordinate
            Else
                Qq(5, J) = Loc_Point(1, J)
            END IF
        END DO
    END DO

    N1 = Panel_Centroid_I(1) - Panel_Centroid_J(1)
    N2 = Panel_Centroid_I(2) - Panel_Centroid_J(2)
    N3 = Panel_Centroid_I(3) - Panel_Centroid_J(3)
    Loc_P(1) = N1*Loc_Coord(1, 1) + N2*Loc_Coord(1, 2) + N3*Loc_Coord(1, 3)      ! (x,y,z)·i
    Loc_P(2) = N1*Loc_Coord(2, 1) + N2*Loc_Coord(2, 2) + N3*Loc_Coord(2, 3)
    Loc_P(3) = N1*Loc_Coord(3, 1) + N2*Loc_Coord(3, 2) + N3*Loc_Coord(3, 3)
    S = 0

    DO I = 1, 4
        C1 = (Qq(I, 1) - Loc_P(1))**2 + Loc_P(3)**2
        C2 = (Qq(I + 1, 1) - Loc_P(1))**2 + Loc_P(3)**2
        H1 = (Qq(I, 1) - Loc_P(1))*(Qq(I, 2) - Loc_P(2))
        H2 = (Qq(I + 1, 1) - Loc_P(1))*(Qq(I + 1, 2) - Loc_P(2))
        IF (Qq(I, 1) == Qq(I + 1, 1)) THEN
            Mm = 0
        Else
            Mm = (Qq(I + 1, 2) - Qq(I, 2))/(Qq(I + 1, 1) - Qq(I, 1))
        END IF
        R1 = SQRT(C1 + (Qq(I, 2) - Loc_P(2))**2)
        R2 = SQRT(C2 + (Qq(I + 1, 2) - Loc_P(2))**2)
        IF (Loc_P(3) == 0) THEN
            S1 = 0
        Else
            S1 = Atan((Mm*C1 - H1)/(Loc_P(3)*R1)) - Atan((Mm*C2 - H2)/(Loc_P(3)*R2))
        END IF
        S = S + S1
    END DO

    Sz = S
    S = 0

    DO I = 1, 4
        L = SQRT((Qq(I + 1, 1) - Qq(I, 1))**2 + (Qq(I + 1, 2) - Qq(I, 2))**2)
        R1 = SQRT((Loc_P(1) - Qq(I, 1))**2 + (Loc_P(2) - Qq(I, 2))**2 + (Loc_P(3) - Qq(I, 3))**2)
        R2 = SQRT((Loc_P(1) - Qq(I + 1, 1))**2 + (Loc_P(2) - Qq(I + 1, 2))**2 + (Loc_P(3) - Qq(I + 1, 3))**2)
        IF (L == 0) THEN
            S1 = 0
        Else
            S1=((Qq(I+1,1)-Qq(I,1))*(Loc_P(2)-Qq(I,2))-(Qq(I+1,2)-Qq(I,2))*(Loc_P(1)-Qq(I,1)))/L*LOG((R1+R2+L)/(R1+R2-L))
        END IF
        S = S + S1
    END DO

    S = S + Loc_P(3)*Sz
    RETURN
END SUBROUTINE Potential

SUBROUTINE Brinv(A, N)
    ! 参  考 ：孙士良《Fortran常用算法程序集》 Gauss-Jordan法实矩阵求逆
    ! 程序编制：左志华
    Implicit None
    Integer(4) :: N
    Real(8) :: A(N, N)                 ! B(N,N),C(N,N)
    Integer(4) :: L                      ! Judge Singular Matrix(奇异矩阵)
    Integer(4) :: I, J, K                ! For Loop
    Real(8) :: T, D
    Real(8), Dimension(N) :: Is, Js

    L = 1

    Do K = 1, N
        D = 0.0
        Do I = K, N
            Do J = K, N
                If (Abs(A(I, J)) .Gt. D) Then
                    D = Abs(A(I, J))
                    Is(K) = I
                    Js(K) = J
                End If
            End Do
        End Do

        If (D + 1.0 .Eq. 1.0) Then
            L = 0
            Write (*, *) 'singer metrix'
            Read (*, *)
            Stop
            return
        End If

        Do J = 1, N
            T = A(K, J)
            A(K, J) = A(Is(K), J)
            A(Is(K), J) = T
        End Do

        Do I = 1, N
            T = A(I, K)
            A(I, K) = A(I, Js(K))
            A(I, Js(K)) = T
        End Do

        A(K, K) = 1/A(K, K)

        Do J = 1, N
            If (J .Ne. K) Then
                A(K, J) = A(K, J)*A(K, K)
            End If
        End Do

        Do I = 1, N
            If (I .Ne. K) Then
                Do J = 1, N
                    If (J .Ne. K) Then
                        A(I, J) = A(I, J) - A(I, K)*A(K, J)
                    End If
                End Do
            End If
        End Do

        Do I = 1, N
            If (I .Ne. K) Then
                A(I, K) = -A(I, K)*A(K, K)
            End If
        End Do

    End Do

    Do K = N, 1, -1
        Do J = 1, N
            T = A(K, J)
            A(K, J) = A(Js(K), J)
            A(Js(K), J) = T
        End Do
        Do I = 1, N
            T = A(I, K)
            A(I, K) = A(I, Is(K))
            A(I, Is(K)) = T
        End Do
    End Do

    Return
End SUBROUTINE Brinv

